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NEEDLE DETECTION, SEGMENTATION AND TRACKING FOR 3D ULTRASOUND- 

GUIDED INTERVENTIONAL PROCEDURES 

Lori Gardi, Donal B. Downey, Aaron Fenster 
Robarts Research Institute 
100 Perth Drive, London, ON, CANADA 



ABSTRACT 

Ultrasound-guided interventional procedures such as breast biopsy and prostate brachytherapy require 
Zt neeXs inserted into the body be tracked in near real-time. In robotxcWed mter^txor^l 
procedures such as robot-aided and ultrasound-guided prostate brachytherapy as well as free-hand 
uWound euided biopsy procedures, needles are inserted free from parallel trajectory constants, i.e., 
Stjlt cTbe positioned with considerable flexibility including oblique trajectories. 
Ho^eT Se insertion results in the needle intersecting the 2D TRUS image and appearing as a 
to bld^uidance. Here, we describe a method for oblique and parallel needle segmenta ion 
S SL be used in 3D ultrasound guided procedures including robot arded prostate 
"achXapy This approach uses near real-time 3D ultrasound (US) to identify segment, tmck and 
d spLf a nSe that is Within the 2D US plane or obliquely intersects the plane of the 2D US image 
Thfs removes the restriction of the needle having to be parallel to the plane of the transducer in current 
WhvAeranv cryotherapy and biopsy procedures. The algorithm applies a grey-level change 
SnteKSe location and orientation of needles from 3D image, Three 2D images 
Sg " e (oblique sagittal, coronal and transverse planes) are £^££S^ 
near real-time. Testing showed that our algorithm's accuracy in segmenting the 3D needle s orientation 
i^ for a chicken tissue phantom, and 0.58° for agar phantoms, over a *15° insertion onentat:on. 
The execution time average is 0.13s on a 1.2 GHz computer. 

1. INTRODUCTION 

This invention has a wide range of applications, such as biopsy of the breast and liver and image 
Sded interventions such as brachytherapy, cryotherapy, as well as other procedures that ^require * 
needle or needles to be introduced into soft tissues and be positioned accurately and Precisely. As an 
example, we describe the use of this invention in robot-aided 3D US-guided prostate brachytherapy. 

Transperineal prostate brachytherapy provides an improved alternative for minimally-invasive 
SeTt of prostate cancer [1] However, pubic arch interference (PAI) with the implant path occurs 
patients with large prostates and/or a small pelvis. These patients cannot be treated with 
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current brachytherapy using parallel needle trajectories guided by a fixed template, because the anterior 
and/or the anterolateral parts of the prostate are blocked by the pubic bone. 

To solve the PAI problems, it would be advantageous to free needle insertions from parallel trajectory 
10 sorveme riup > positioned with considerable flexibility including oblique 

SSS V 3] t e oX-Ta 1- patients with PAI to be treated with brachytherapy 

Ssf f d Jgoing lengthy hormonal downsizing therapy In addiflon, ci^W 
size during the implantation may require re-optimization of the dose plan. Thus, combming mree 
dtoenZal (3D) transrectal ultrasound (TRUS) imaging, dosimetry and obhque needle inserts 
Soriefcan provide the tools needed for dynamic re-optimization of the dose plan dunng seed 

by allowing dynamic adjustments of the needle position to target potential 

"cold spots". 

To achieve this goal, we developed a 3D TRUS-guided robot-aided prostate brachytherapy system 
Ml to ttL^S a robot is used as a movable needle guide containing one-hole needle guide 

and orientation of the needle can be changed as the robot moves. By unrfymg 
he robo ultrasound (US) transducer and the 3D TRUS image coordinate systems, the position of the 
temptte hole can be accurately related to the 3D TRUS image coordmate system all — 
and consistent insertion of the needle via the hole into the targeted position in the prostate along 
various trajectories including oblique ones. 

However, oblique insertion will result in the needle intersecting the two-dimensional TO.TRUStajB 
Pl Je te the needle will only appear as a dot in this image (see Figure 2). Some investigators have 
a"d automatic needle sedation metho ds to locate needles for biopsy and ^thenmy [5,6]. 
However, these methods require that the needle be completely contained in the 2D US image [6]. 

In this paper we describe a near real-time method, for automatic identification, segmentation and 
Ltmg If neX during oblique insertion even if the needle exits the 2D US nnage plane This 
SSd can aTso be used for automatic identification, segmentation and tracking of needles if they 
completely contained in the 2D US image plane. 

2. EXAMPLE APPLICATION 
2.1. 3D TRUS-guided robotic-assisted system (see Figure 1) 

The 3D TRUS imaging system used to test the needle segmentation method co * sist *^ 
12 GHz processor for 3D image acquisition, reconstruction and display, a video frame grabber for 
^^"ition, a mov?r conLler module (MCM), which controls the rotation of the trans 
recS uEJd (TRUS) transducer-mover assembly via the serial port of the computer. To produce 
the MCM and motor assembly causes the transducer to rotate about its long axis 
over about 100°while 2D TRUS images are digitized at intervals by the frame grabber [7]. These 
acquired images are reconstructed into a 3D image and are available for viewing in real-time as the 2D 
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images are being acquired. The 3D TRUS image is viewed using 3D visualization software, which 
includes multi-planar 3D display as well as an extensive set of measurement tools [10-12]. 

The robot svstem (see Fig. 1) includes a robotic arm assembly with 6 degrees-of-freedom. The 
positioning software can control the robot arm to move in terms of world/tool coordinate systems The 
world coordinate system is fixed to the ground and the tool coordinate system is fixed to the robot 
« The robot control software has been integrated together with the 3D US visualization software 
T4 7] so that we can perform precise image-based robot path planning and visualize multiple potential 
trajectories in the 3D TRUS image coordinate system using a graphical user interface. Figure 1 shows 
the whole system. 

2.2. Needle segmentation algorithm 

To perform near real-time needle segmentation, capture of two 3D US images is required. Note that if 
the needle is contained in a 2D image, then two 2D US images are required, but the procedure , is 
unchanged. The first (pre-scan) 3D image is obtained by scanning the prostate (tissue) before he 
needle is inserted, and the second (post-scan) is acquired by scanning only the region conning Jie 
needle during needle insertion. The second 3D image is compared against the first and the needle 
position within the post-scan 3D image, including entry point and needle tip location, is determined 
using the approach described below and the flow chart shown in Fig. 5. 

The needle segmentation method can be used for both linear and rotational 3D US scanning approaches 
171 In addition, our segmentation method can be applied equally well to 3D US images reconstructed 
ushig the linear scanning geometry, but acquired using rotational 3D scaiming geometry such as that 
used in prostate imaging. Figure 4(a) shows 3D US images acquired and reconstructed usuig tiie 
rotational geometry approach, but Fig. 4(b) was acquired in the rotational approach and recounted 
using a linear scanning 3D geometry [7]. In Fig. 4(a) needles are stretched horizontally in the 3D 
scanning direction at the center of the image but are more angulated away from the center. In Fig. 4(b) 
needles are spread horizontally everywhere in the image, which greatly simplifies the identification and 
processing image voxels that contain the needle. 

2. 2.1. Grey-level threshold value (see Figure 6) 

A grey-level threshold value is used to reduce the number of voxels to be analyzed and to obtain 
candidate needle voxels. To determine the threshold value, the maximum grey-level value GL max in the 
post-scan image is first found. The threshold value is determined by: 

Threshold value = a x GL max ^' 

where 0 < a < 1. In our testing, we used a = 0.5. 

2.2.2. Generating 3D difference map(see Figure 7) 

Voxels with a grey-level change that exceed the threshold value between pre-scan and post-scan are 
assumed to be part of the needle. In our system, a candidate voxel (CV) containing the needle is defined 
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as a voxel in the difference 3D image between the pre- and post-scans that exceeds a pre-defined 
threshold: 

GLC(ij,Jfc) = postGLOV,*) - preGLftM) < 2) 

where P reGL(U*) and postGLft/,*) are the grey-level values at voxel location QJ.k) inboth the pre- 
scan and the post-scan respectively, and GLC(iJ,fe) is the grey level change. A collection of these 
candidate voxels is referred to as the Difference Map (DM), i.e., 

(31 

(f B j w ye3DDM K) 

where GLC(w^) > Threshold Value; m = 1,2,..., n; and « is the number of points included in the 
3D difference map (see Fig. 9(a)). 

2.2.3. Selecting needle points (See Figure 8) lfc< . 

In our system! we have the advantage that the 3D image is reconstructed on demand, 

always have Access to the original acquired image data. Since the unage of the needle is spread 

horizontally as shown in Fig. 3 (in the 3D scanning direction), we can remove those voxels from our 

DM that are not connected to other voxels in this direction as follows: 

(W«.U£3D DM 
when \jGLC(i m J m ,k m ±s) < threshold value 

where s = 1 2 P/2: and P is the number of voxels surrounding voxel (r ra , j m , k m ) in ^-direction. We 
cTseP = 4 in our algorithm. Figure 9(b) is an example of the DM after the spurious voxels have been 
rejected. 

2.2.4. Finding the needle trajectory (Block 440 in Figure 5 and Figure 10) dl 
At this stage needle candidates are available. There are a number of means for identifying the needle 
Rectory from these candidates. For example, the well-known Hough Transform technique [13] or 
linear regression may be used, as shown in Blocks 444 and 446 of Figure 10. 

Hough Transform approach: We have developed a fast implementation of the Hough Transform [14], 
wS we called the Real-Time Hough Transform. Due to the robustness to the addition of extraneous 
noise the Hough Transform is one of the most powerful line detection techniques nowadays and has 
been widely used in different areas. Unfortunately, its high computation needs often prevents it from 
being applied in real time applications without the help of specially designed hardware Our fast 
implementation of the Hough Transform is based on coarse-fine search and the determination jo the 
optimal image resolution. Comparing to conventional techniques, our approach decreases the tune for 
needle segmentation about two orders of magnitude. Experiments with agar phantoms and patient 
breast biopsy US image sequences showed that our approach can segment the biopsy needle m real- 
time (i.e. less than 33ms) on an affordable PC computer without the help of specially designed 
hardware with the angular rms error about and the position rms error about 0.5 mm. 
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Linear regression approach: An alternative approach is to use linear regression to find the trajectory 
oTZZZte Tsee pfg^re 10). The voxels in the difference map (DM) can be fit to a kne using ; linear 
^nce the equation of the line (needle) is obtained, the difference map is further 
P r?ces^ea to remove voxels that are far away (> 5 voxels) from the line and linear regression is 
perfoLd again with the filtered data. This procedure is repeated until all the voxels in ^the ; differ^ce 
^arewitLapres^ Figur^e 9(c) is an example of the DM before 

me final regression calculation. Once the equation of the line/needle is found, the trajectory of the 
nTedte is known. The needle entry point is the intersection of needle trajectory and the known entry 
plane. The needle tip is the farthest needle voxel along the needle trajectory. 

225 Extracting of the 2D planes containing the needle 

To extract any plane containing the needle, we use the segmented needle entry point, needle tip pom 
and we select a third point within the 3D image to define a specific plane. In this way, we can jextrac 
^display oblique saggital, coronal and transverse views with the needle trajectory highlighted (as 
Ascribed in the flowed of Fig. 3). Figure 11 shows the result of this complete proced ure ^obtained 
during a patient's prostate cryotherapy procedure, demonstrating that the needle can be tracked as it is 
being inserted and orthogonal views can be displayed for the user during the insertion procedure. 



2.3 Evaluation 



2 31 Experimental apparatus 

The accuracy and variability of our needle segmentation and tracking technique was tested using images 
acquired by scanning phantoms. The robot was used to insert a needle at known angles, including 
oblique trajectories with respect to the TRUS image plane (Fig. 1). 

The needle used in these experiments was a typical 18-gauge (1.2mm diameter) prostate 
brachytherapy needle. The US tissue-mimicking phantoms were made of agar, using ^Jf? 
by Ricky et al [8], and chicken breast tissues [9]. TRUS images were obtained usmg ; ar t 8558/S side- 
tog Unear array transducer with a central frequency of 7.5MHz, attached to a B-K Medical 2102 
Sc US machine (B-K, Denmark). The 3D TRUS system consisted of a Pentium III personal 
cZputer (PC) equipped with a Matrox Meteor II video fiame grabber for 30Hz video image 
acquisition. 

2 3 2 Algorithm execution time u 
Execution time depends on the 3D scanning angular interval and the extent of the region to be 
investigated. To evaluate the algorithm's execution time, first we performed the pre-msertion scan, 
and then inserted the needle. After needle insertion, the phantom was scanned again, and the needle 
was segmented. A software timer was used to measure the time duration of the segmentation 
algorithm. 
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ToLiZZtZ of the algorithm, we used the robot to guide the needle insertion into the phantom 
at known angles (Fig. 1). The angulation accuracy of the robot was evaluated to be 0.12 ± 0.07 [4]. 

First we used the robot to guide the needle insertion along a trajectory parallel to the transducer, 
2 the 0* oLtatiofsince the needle could be verified by observing the needle m the real- 
time 2D TRUS image, we assumed that this trajectory was correct. Thus, our obhque trajectory 
Curacy measurement; could be made with respect to the 0° trajectory. The portico, of ^needl Up 
and ent^ point were then found for 0° using the algorithm described above. Then, the robot was used 
to inserTthe needle at different angles (.5>U».S-) with respect to the 0° 

each insertion, the positions of needle tip and entry point were found. The corresponding segmented 
needle vectors through these entry and needle tip points were determined by: 

cos0 " 8 " |a||b| 

whereAis the segmented needle vector for 0° insertion; Bis the segmented needle vector for the 
insertion at any olher angle; 0 alg is the angle derived from the segmentation algonthm. The accuracy of 
the algorithm was evaluated by comparing 0 alg with the robot orientation angle 0 rob , i.e., 

The accuracy test was repeated with a chicken tissue phantom, and the accuracy was also determined 
^ ~) -d (5). For the agar phantoms, we performed five groups of tests to evaluate the 
algorilhm execution time and accuracy. Each group consisted of seven insertions ,i.e., insertion at 
0" and ± 5',±10',±15- . The mean error as a function of insertion angle was calculated, thus: 

em = i|(0. te ),-<0.*)i|/3- (6) 

3. RESULTS AND CONCLUSION 

Table 1 shows the evaluation results. In the chicken tissue phantom, the average execution time was 
0.?3 im seconds, and the average angulation error was 0.54« f 0.16°. £ ^£ 

execution time was 0.12 ±0.01 seconds, and the average angulation error was 0.58 ± 0.36 . Th > results 
shown in Table 1 also demonstrate that the insertion error does not significantly depend on insertion 
angle. 

In the 3D US images, needle voxels generally have high grey level values However due to ^pecular 
reflection, some background structures may also appear to have high grey level values 
the difficulty in automatic needle segmentation in an US image using grey level information directly 

6 

Needle segmentation patent disclosure 



K„.i»DTn <>nm the IhW Ima ge Database on OZHb/gOOS 



1/12/04 



rs 61 Because US images suffer ftom low contrast, signal loss due to shadowing, refaction and 
^aZarrffte'aW.eve. change detection techni q ue appears to * addrfon, 

since the needle is segmented from a difference map, complex backgrounds are removed. 

,» conclusion, we developed a ^ey-ieve. change detection lechniqu e and tesKd tatajbUjr £™ 

SSTJS^ SSI approach hafalso been tested during severa. prostate cryotherapy 
orocedures with excellent results. 
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Figure 1. 3D TRUS-guided and robot-assisted prostate brachytherapy system. 
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2D Images 




1 1£ Transducer 



Figure 2. Schematic diagram showing the ultrasound transducer mounted on the Rotational 
moused to produce 3D TRUS images of the prostate. The diagram shows tha when the 
needle is inserted at an oblique angle relative to the transducer's axis, multiple 2D US images are 
needed to track its insertion trajectory. 
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Insert needle 
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NEEDLE 



Calculate Needle 

Trajectory 
(from needle seg.) 



Calculate needle 
entry point and 
needle tip locations 
within raw volume 



Reconstruct the 


needle tip and entry 
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Extract the plane 
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reconstructed image 




f 


Display the 


extracted plane 






Repeat for other 


planes. 





Figure 3. How chart showing the complete needle identification, segmentation method and 
display. The step labeled "SEGMENT NEEDLE" is detailed in Figure 5. 
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410 



Calculate GLC 
threshold 




I 



420 



Generate Difference 
Map (DM) from raw 
3DUS image 




440 



Perform regression 
analysis (eg linear) 
on DM 




460 



Filter DM 




NO 




Figure 5. Details of Block 400 "SEGMENT NEEDLE" in Fig. 2, i.e. the steps in segmenting 
the needle using acquired 2D or 3D US images. 
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Init: GLmax = 0 



I 



Select next voxel 
from 3DUS image 



No 




Get grey-level value: 
GLvalue 



No 




Yes 



Threshold = 
GLmax * a 



GLmax = GLvalue 



Figure 6. Flow chart for calculating the GLC threshold in Block 410 in Fig 5 Analyze all pixels 
in the 3DUS post-scan and find the maximum grey-level (GLmax) value. I is assumed that he 
needle data will have the brightest signal and will generate the highest voxel values in the 3DUS 
image. The threshold value selected must be less than the maximum grey-level value in the post- 
scan image. We calculate the threshold value as follows: 

Threshold value = GLm * a; where 0 < a < 1 



where optimum value of a is found by experimentation. 
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Select next voxel in 
the 3DUS image 



Find grey-level 
difference (GLdiff) 
(post-pre) 



No 




Yes 



Add voxel position 
to difference map 
(DM) 



Figure 7. Details of Block 420 in Fig. 5. Flow chart for generate the difference map (DM) from 
the 3D US image. A difference map (DM) is an array of 3D voxel locations (x,y,z) within the 
3DUS image where the grey-level difference between pre and post scan exceeds the CjLC 
threshold value. If 2D images are used, then DM is a 2D map. 
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Firnire 8 Details of Block 460 in Figure 4, in which the Difference Map (DM) is filtered. After 
dSLg trajectory of the needle, we throw away points in the DM that are far away (by n 
pS ftom the equation of the line defining the needle. We then perform a least squared fit 
(LSF) on the remaining points in Block 480. 
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Fieure 9 Results of the 3D difference map (DM) and processing procedures, (a) 3D difference 
S (DM) I thout any processing; (b) the 3D difference map after spurious s.gna is have been 
rejJctedrand (c) final 3D difference map (DM) after spurious voxels have been rejected. 



Figures for needle segmentation patent disclosure 



9 



1/12/04 




442 



1 f 

Select ne 
the 


xt point in 
DM 


i 


r 


Project volume 
along z~axis and find 
trajectory 1 




f 


Project volume 
along y-axis and find 
trajectory 2 




444 




446 




No 




442 



Given X, use results 
of trajectory 1 to get 
Y and trajectory 2 to 
getZ 




Figure 10. Details of Block 480 in Figure 4, in which a least squares i ffl . is performed I on fte 
pixels in the DM. The collection of 3D points (x,y,z) referred to as the difference map (DM) 
contains the voxel positions in the 3DUS image that potentially belong to the set of points 
deSg the needle P To determine the trajectory of the needle, we first apply a 2D least: squares 
fit (LSF) on the x,y component of the voxels in DM. Then, we apply the same 2D LSF to fte x z 
component of the voxels in DM. The result, are two perpendicular projections of the needle in 
2D space, which when combined, give us the trajectory of the needle m 3D space. 
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Figure 11. Example of Oblique needle segmentation needle tracking obtained during a 
cryosurger; procedure on a patient's prostate, (a) Oblique sagital view; (b) obhque coronal 
plane; (c) transverse view with needle projected 
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Table 1 Test results for chicken tissue and agar phantoms. 1. Test results for chicken tissue 
phantom; 2. Test results for agar phantoms. 



Angle 


-15* 


-10* 


-5' 


5 - 




15- 




Time (S) 


0.13 


0.11 


0.12 


0.12 


0.12 


0.14 


1 


Accuracy 


0.50° 


0.51° 


0.43° 


0.37° 


0.74° 


0.74* 




Time (S) 


0.12 


0.12 


0.12 


o.u 


0.12 


0.13 


2 


Accuracy 


0.30* 


| 0.71° 


0.48° 


0.68° 


0.42° 


0.86° 
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